rm(list = ls())
library("ROCR")
library("magrittr")
library("xtable")
load("results/pr_fcast_agg.rda")
load("results/pr_test6.rda")

####################################################################################
############# Table 4
####################################################################################
pr_test6 <- dplyr::filter(pr_test6, !is.na(y))

# p >= cutoff -> 1
t_select <- function(type=NULL, target=NULL, pred, obs) {
        if (is.null(type)) stop("Specify type: ?performance")
        pred_obj <- prediction(pred, obs)
        perf_obj <- performance(pred_obj, type, x.measure = "cutoff")
        stats   <- perf_obj@y.values[[1]]
        cutoffs <- perf_obj@x.values[[1]]  # >= is 1
  # Find closest value to desired target; print message if not = target
  idx <- which.min(abs(stats - target))
  if (!isTRUE(all.equal(stats[idx], target))) {
    message("Closest value to target is ", stats[idx])
  }
  cutoffs[idx]
}
# For a given cutoff, find performance
t_fit <- function(cutoff, want, pred, obs) {
    # Convert to factor so we don't lose 0 entries in confusion matrix
    obs <- factor(as.numeric(obs), levels=c(0, 1))
    pred <- factor(as.numeric(pred >= cutoff), levels=c(0, 1))
    cmat <- table(obs, pred)
    if (want=="tpr" | want=="rec") {
      res <- cmat[2, 2] / sum(cmat[2, ])
    } else if (want=="fpr") {
      res <- cmat[1, 2] / sum(cmat[1, ])
    } else if (want=="prec") {
      res <- cmat[2, 2] / sum(cmat[, 2])
    } else if (want=="fnr") {
      res <- cmat[2, 1] / sum(cmat[2, ])
    } else {
      stop("Unrecognized want '", want, "'")
    }
    res
  }
  
  
t1 <- t_select("tpr", 0.5, pr_test6$ebma, pr_test6$y)
message("Recall: ", round(t_fit(t1, "rec", pr_test6$ebma, pr_test6$y), 2))
message("Precision: 1 in ", round(1/t_fit(t1, "prec", pr_test6$ebma, pr_test6$y)))
message("FPR: ", round(t_fit(t1, "fpr", pr_test6$ebma, pr_test6$y), 2))
  
# FNR is 1 - recall

tab3 <- pr_fcast_agg %>%
        dplyr::filter(ebma >= t1) %>%
       # mutate(country = prettyc(country)) %>%
        dplyr::select(gid, ebma) %>%
        xtable(digits=3, label="top-fcast") 
  
# Higher recall, broader forecasts
t1 <- t_select("tpr", 0.76, pr_test6$ebma, pr_test6$y)
t_fit(t1, "prec", pr_test6$ebma, pr_test6$y)
t_fit(t1, "fpr", pr_test6$ebma, pr_test6$y)

pr_fcast_agg %>%
       dplyr::filter(ebma >= t1) %>%
        dplyr::select(gid, ebma) %>%
        xtable(digits=3)
save(pr_fcast_agg, file = "results/pr_fcast_agg.RData")

load("data/pr_fcast_agg.RData")

library(cshapes)
pr_fcast_agg_grid <- pr_fcast_agg %>% 
            dplyr::slice_max(order_by = ebma, n = 15) %>% 
          dplyr::pull(gid)
gid <- pr_fcast_agg_grid
# produce a map for myanmar only
myanmar_map <- cshp(date = as.Date("2010-12-31"), useGW = T)
myanmar_map <- myanmar_map[which(myanmar_map@data$GWCODE==775),]

myanmar_map2 <- fortify(myanmar_map, region = "COWCODE")

load("data/myanmar_priogrid_cell.RData")
# select four gids
cooridates <- myanmar_priogrid_cell %>% 
          filter(gid %in% pr_fcast_agg_grid) %>%
          dplyr::select(gid, xcoord, ycoord)

cooridates_df <- cooridates %>%
  filter(gid == 158947 | gid ==159667 | gid == 158946 | gid == 159666 | gid == 162556)

library(ggrepel)

ggplot() + 
 geom_polygon(data = myanmar_map2, 
               aes(x = long, y = lat, group = group) , 
               fill = "gray") + coord_fixed() +
  geom_point(data = cooridates, aes(x = xcoord, y = ycoord), 
             shape = 15, size = 6, color = "blue") +
  geom_text_repel(data = cooridates,
                  aes(xcoord, ycoord,label = gid, group = gid),
                  color = 'blue',
                  max.overlaps = 15,
                  size  = 6,
                  box.padding = 0.4, point.padding = 0.5) +
  scale_colour_distiller(type = "seq", palette = 1, na.value = "gray", guide = "colourbar") +
  scale_fill_manual(values = c("white","gray"), 
                    name= "", guide = guide_legend(reverse = TRUE)) + 
  labs(title = "",
       subtitle = "",
       y = "纬度（北纬）", x = "经度（东经）")+theme_bw()+
  theme(axis.ticks.length = unit(0,"cm"),
        panel.spacing = unit(0, "lines"), 
        plot.title = element_text(hjust = 0.5),
        text = element_text(size=15, family = 'KaiTi'),
        plot.margin = unit(c(0,  0, 0, 0), "lines")) 
ggsave("figs/table_4a.png", width = 6, height = 8)

agg_grid <- pr_fcast_agg %>% 
      dplyr::slice_max(order_by = ebma, n = 15) %>% 
      dplyr::select(gid,  ebma)

agg_grid <- left_join(agg_grid, cooridates, by = "gid")
library(xtable)  
agg_grid$geometry <- NULL
agg_grid <- as.data.frame(agg_grid)
xtable(agg_grid,digits = 6) %>% 
    print(file = "figs/table4.html",type = "html")
####################################################################################



